Diffusion Monte Carlo with lattice regularization 



Michele Casula, 1,2 Claudia Filippi, 3 and Sandro Sorella 1,2 
1 International School for Advanced Studies (SISSA) Via Beirut 2,4 34014 Trieste, Italy 
2 INFM Democritos National Simulation Center, Trieste, Italy 
3 Universiteit Leiden, Instituut-Lorentz for Theoretical Physics, NL-2300 RA Leiden, The Netherlands 

(Dated: February 2, 2008) 

We introduce an efficient lattice regularization scheme for quantum Monte Carlo calculations of 
realistic electronic systems. The kinetic term is discretized by a finite difference Laplacian with 
two mesh sizes, a and a', where a'/a is an irrational number so that the electronic coordinates 
are not defined on a particular lattice but on the continuous configuration space. The regularized 
Hamiltonian goes to the continuous limit for a — + and provides several advantages. In particular, it 
allows the inclusion of non-local potentials in a consistent variational scheme, substantially improving 
the accuracy upon previous non- variational approaches. 
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In the last few decades, enormous progress has been 
made in computing the physical properties of many- 
electron systems with numerical methods based on first 
principle quantum mechanics. In particular, quantum 
Monte Carlo (QMC) techniques Q have proven very suc- 
cessful mainly because they allow the explicit inclusion 
of electronic correlation in a many-body wave function 
(WF). Moreover, projection QMC methods such as dif- 
fusion Monte Carlo (DMC) can improve upon a given 
guiding WF $g by stochastically projecting a state ^fn 
much closer to the exact ground state (GS) of the Hamil- 
tonian H. ^fn is obtained within the fixed- node ap- 
proximation (FN A), which yields the lowest solution of 
the Schrodinger equation with the same nodes as f g- If 
^Pg is appropriately optimized, the method usually pro- 
vides a good upper bound Efn for the GS energy Eos'- 
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In general, in a QMC calculation, one accesses a mixed- 
average estimate of the total energy: 
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which coincides with Efn since <J fn satisfies the equa- 
tion fn — Efn^fn within the nodal pockets of ^q- 
While the QMC methods can be extended to large sys- 
tems containing many electrons, the computational ef- 
fort increases dramatically with the atomic number Z. 
The most common way to overcome this difficulty is to 
replace the core electrons by pseudopotentials, an ap- 
proximation which is usually rather good as the core is 
chemically inert. Since the pseudopotentials are in most 
cases non local, H contains a non-local potential V p ', 
i.e. (x\V p \x') / even when \x — x'\ ^ 0, where x de- 
notes an all electron configuration with positions {r^} 
and spins {(Ji\. Consequently, the standard DMC ap- 
proach cannot be applied and the so called "locality ap- 



proximation" (LA) is employed 0, EI El which ap- 
proximates the non-local potential V p with a local one, 
V LA (x) = (x\V p \$>g) I '{x\^g}- A major drawback of 
this approach is that Efn (Eq. is no longer available 
since it does not coincide with Em a (Eq. EJ) accessible in 
a DMC calculation: ^fn is now the best solution with 
the same nodes as for the approximate Hamiltonian 
with potential V LA and not for the Hamiltonian H. Even 
though Em a equals Eqs if is exact, it is otherwise 
not bounded by Eqs and no rigorous information about 
the quality of ^fn is given by Em a- For instance, it 
is not possible to accertain the quality of the nodes of 
since a lower Em a may correspond to a worse ^fn, 
namely with higher Efn- Empirically, the method ap- 
pears to work , but its drawback has limited the impact 
of this technique to a wider range of applications. 

In this Letter, we present a lattice regularization of 
the many-electron Hamiltonian which removes the above 
difficulties when using non-local potentials within the 
FNA. We demonstrate the utility of our lattice regular- 
ized DMC approach in cases where the LA yields inac- 
curate results. 

Regularization of the Hamiltonian. We consider the 
Hamiltonian in atomic units: 
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and, for the moment, assume a local potential, V(x) = 

J2 3 - z j/\n - Rj\ + J2i< 3 Vl^ _ fifl) wnere R an d Vin- 
dicate the ionic and electron positions, respectively. 

The Laplacian is approximated by a finite difference 
form 
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where A a ' p is an appropriate Hermitian lattice operator 
defined by a mesh size a, a constant r\, and an arbitrary 
function p(r,): 

K" P f( x i,Vi, z i) = v/a 2 [p(x l + a/2)(f(x l + a) - f{xi)) 
-I- p(xi - a/2)(f(xi - a) - /(a*))] + Xi <-> yi <-> z { . 



2 



For p = i] = 1, A°' p coincides with the usual discretized 
form of the Laplacian on a lattice with mesh size a. By 
combining two such operators in Eq.0|with mesh sizes a 
and a', and a' /a an irrational number, the electron co- 
ordinates {n} are still defined in the continuous space 
because the two meshes are incommensurate and the dif- 
fusion process described by this discretized Laplacian can 
cover all the configuration space. Moreover, in Eq.QJ the 
function p can be arbitrarily spatial dependent as long 
as < p < 1 to ensure that p and 1 — p are both posi- 
tive. In order to optimize the efficiency of the diffusion 
process, we have carefully chosen the constant a' /a and 
the function p as: 



p(r) = 1/(1 + Z 2 /4\f~R\ 2 ). 
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where R and Z are the position and the atomic number of 
the nucleus closest to the electron in r. Therefore, if the 
electron is very close to a nucleus, the smaller mesh a is 
used, while, in the valence region where \r — R\ » 1/Z, 
the electron can make steps of larger amplitude a' , thus 
reducing the QMC correlation time. 

If 1] = 1 + 0(a 2 ), the resulting finite difference Lapla- 
cian A" coincides with the continuous A up to 0(a 2 ). 
To further improve the accuracy of the approximation 
and work with reasonably large value of a, we regularize 
also the potential V — > V a so that our final regularized 
Hamiltonian, H a = -1/2^ A? + V a , fulfills the follow- 
ing three conditions: i) H a — > H for a — > 0; ii) for the 
chosen guiding WF for any a and any x, the local en- 
ergy cl(x, [v&g]) = (xlH^a) / (x\^g) of the continuous 
Hamiltonian H is equal to e1(x, [v&g]) corresponding to 
the Hamiltonian H a ] iii) the discretized kinetic energy is 
equal to the continuous one calculated on the state ^>g, 
i-e- (-1/2 Ei A<V g = (-1/2E,A,) $ ,. While the con- 
dition (iii) determines the constant rj, the condition (ii) 
leads to the following solution for V a : 



V a (x) = V(x) + 
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Notice that the condition (ii) yields another important 
property for H a : if is an eigenstate of H, it is also an 
eigenstate of H a for any a, as can be derived by simple 
inspection using that A a is Hermitian. Thus, by improv- 
ing ^g, a better a — > convergence is also expected. 

Non-local pseudopotential. The above regularization 
also applies when we include in the potential V a non- 
local potential of the form V p = E, vP (^i) where 
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The sum over j is here restricted to the pseudoatoms, 
Pi is the Z-angular-momentum projection centered on the 
pseudoatom at Rj , and V max the corresponding maximum 
angular momentum considered. The functions Vj and i>j 
are radial, and vanishes outside a core radius rj_. 

In a QMC calculation with pseudopotentials 0, Q , the 
angular integration to evaluate the projection Pi is per- 
formed numerically on a regular polyhedron defined by 
Ny vertices and by employing a quadrature rule for the 
integration on a sphere. Therefore, the non-local pseu- 
dopotential v p acts on a configuration x by means of 
a finite number of matrix elements equal to NvN core , 
where N core is the number of electrons within the core 
radius of a pseudoatom. In the following, we will assume 
that, in the continuous Hamiltonian H, the pseudopoten- 
tial is by definition discretized with a finite number Ny 
of vertices. Then, the regularization H — » H a follows by 
simple substitution of V(x) — > V(x) + V p in Eq. H3 

Lattice regularized diffusion Monte Carlo (LRDMC). 
Although H a is an Hamiltonian defined on a continuous 
space, all techniques valid on a lattice can be straight- 
forwardly applied here since H a acts on a configuration 
exactly as a lattice Hamiltonian, namely: 



|ff a |*G} = $^V<z'|*G> 
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where, for a given x, the number of matrix elements 
if" x , are finite and can be computed even in the pres- 
ence of non-local pseudopotentials. In particular, we can 
resort to the same scheme used in the efficient lattice 
Green function Monte Carlo algorithm 0, 0, E|- The 
resulting algorithm, valid for a continuous regularized 
Hamiltonian, is the LRDMC. The corresponding Green 
function matrix elements are G x >. x = ^ q(x' ){K5 x i ^ x — 
H%, ^/^cix), and, provided they are all non-negative, 
the positive distribution ^/g(x)^gs(x) is statistically 
sampled. Note that, since the spectrum of H a is not 
bounded from above, we need to take the limit A — > oo, 
which can be handled with no loss of efficiency as de- 
scribed in Ref. 0. The LRDMC algorithm is very sim- 
ple: given a walker with configuration x and weight w, 
a new configuration x' ^ x is obtained with probability 
Px'.x = G x ' tX /N x , where N x = J2x>(=£x) G x',x is the nor- 
malization. The walker weight is then changed by a fac- 
tor w — > w exp(-T x e L {x, [*g])) where t x = -log(r)/N x 
is a diffusion time determined by a single random num- 
ber < r < 1. Then, the usual branching scheme is used 
to control the fluctuations of the walker weights. 

Within this approach, the time-step error of the stan- 
dard DMC method is not present, but is replaced by 
the error of the Laplacian discretization. By performing 
LRDMC simulations of H a for different et's, we can ap- 
proach the GS properties of a continuous system in the 
limit a — * 0. The method becomes less efficient as a — > 
since the computational time increases like 1/a 2 , simi- 
larly to the DMC efficiency which is limited by 1/Ar, 



3 



CD 



X 



-37.832 



-37.838 



ra -37.844 



-37.85 



DMC 
LRDMC 



0.01 0.02 

time step (Hartree 



0.03 

1> 



0.04 



FIG. 1: FN energies for the all-electron carbon atom com- 
puted within DMC and LRDMC. The lattice spacing a has 
been mapped to the time-step r using the relation a — a/t. 

where At is the imaginary time step for the diffusion 
process. 

Since the Green function G x >, x can be made strictly 
positive only for bosons, we have to introduce here the 
analogous of the FNA on a lattice 0, 0- by modify- 
ing few of the matrix elements of the regularized Hamil- 
tonian H a . For each configuration x, the matrix el- 
ements H£, x which yield G X > :X < are set to zero 
and included in the so called sign-flip term, V s /(x) — 
12x'jtx'^G( x ') H x',x/^G(x) > 0, which is then added to 
the diagonal element Hi 1 x Q. The resulting effective 
Hamiltonian H cS has the same local energy as H a , it is 
non frustrated, and its GS WF has the same signs as the 
guiding function ^>c(x). The GS energy of H eS can be 
efficiently computed with the mixed average Ema and, 
for a local Hamiltonian H, we recover the standard DMC 
result Em a — Ep n in the limit a — > as shown in Fig. ^ 
Moreover, when H contains non-local pseudopotentials, 
Em a for the GS of H eS is certainly lower than the ex- 
pectation value Eq of the Hamiltonian H a on "F^, and 
represents an upper bound of the variational energy Ep jy 
with the Hamiltonian H a in Eq. ^ as was shown for the 
lattice FN approach 0. This important upper bound 
property has been so far only obtained within the lattice 
regularization proposed here, whereas, in the DMC ap- 
proach, it is not true in general that Eq > Em a > Efn- 

Computation of Epn- Another advantage of our lat- 
tice regularization is that we can also compute directly 
and accurately Ep jy for the Hamiltonian H by introduc- 
ing a more general effective Hamiltonian characterized by 
the parameters 7 > and < a < 1 10]: 
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H a x , x + (1 + i)V. f {x) + a(l + 7 )V s p / (a;) (9) 
-jH^ x if * G {x')H a xltX /^ G (x) > 
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where x' ± x and V?{x) = *gV)V£J*g{x) < 

0. This Hamiltonian also satisfies G x / yX > and reduces 
to the standard lattice FN one for a = 7 = 0. 

The parameter 7 allows one to evaluate Epjq in the 



presence of the non-local contribution of V p as 

dE M A{l) 



E FN = E M a(j) - (7 + 1)- 
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where we used that H — H cS — (7+ l)d 7 H eS and applied 
the Hellmann-Feynman theorem to the last term. To es- 
timate Epn, we compute the derivative with respect to 
7 in an approximate but variational way, the best varia- 
tional estimate being for 7 = [Ioj |: 

E FN < E M a{0) - [E MA {l) - E MA {0)}h, (11) 

where the equality sign holds in the limit of small 7. No- 
tice that, for a = 1 and 7 = 0, the LA is recovered, i.e. 
all the off-diagonal matrix elements of V p are included 
in the diagonal ones of H cS . Even though it is not garan- 
teed that Em a is variational for a > 0, with the present 
scheme we can evaluate Epjy and improve upon the lo- 
cality approximation by optimizing a. Indeed, the vari- 
ational energy Ep^ corresponding to a given value of a 
can be estimated efficiently using correlated sampling to 
evaluate Em a for 7 = and < 7 < 1/a — 1 as in 
Eq. ^2 Remarkably, this approach allows us to work 
close to the locality approximation, even with a = 0.95 
and 7 = 0.05, with statistical errors for Ep^ similar to 
the ones for Ema- 

Results and perspectives. We have first tested the per- 
formance of the LRDMC approach on the silicon pseu- 
doatom using three Hartree- Fock pseudopotentials which 
differ in the construction, functional form and core ra- 
dius. For each pseudopotential, we employ three WF's 
with the same determinantal component and, conse- 
quently, the same nodes, but with different Jastrow fac- 
tors. We use no Jastrow factor, a two-body, and a so- 
phisticated three-body Jastrow factor . 

As shown in Fig. |3 the energy estimate Ema com- 
puted within DMC changes significantly with the guid- 
ing WF ^>g, and differs from the variational expectation 
value, Epn, which we compute with the LRDMC scheme 
for a = 0, 0.5 and 0.9. For all cases, the statistical un- 
certainty does not allow us to discriminate between the 
LRDMC energies obtained for a = 0.5 and a — 0.9, and 
a shallow minimum seems to lie between these two val- 
ues. The optimal Ep^ is almost free of the localization 
error and the weakest de pen dence on ^>g is obtained for 
Lester's pseudopotential [1J| which has the smallest core 
radius in the non-local component. Interestingly, since 
Efn for oi ~ 1 is very close to the minimum, the LA 
seems to yield good WF's. However, the value of a should 
be optimized for each case since there is no general rea- 
son why a = 1 should give the best energy. In principle, 
in the case of systems with lower symmetry, the value of 
a should be closer to due to the non conservation of 
the total angular momentum. 

A stringent test case for our LRDMC algorithm is 
the scandium atom: the LA for transition metals yields 



4 




CD 
CD 

r 

I 

>. 
a> 

cd 
c 

CD 



0.02 0.04 0.06 0.08 

VMC energy - FN energy (Hartree) 

FIG. 2: FN energies of the silicon pseudoatom computed 
within DMC and LRDMC. For different pseudopotentials 
(Soft Dolg's 0| and Lester's 0|), we use as guiding 

WF's a Hartree-Fock determinant with no Jastrow, a two- 
body and a three-body Jastrow factor. A more accurate 
guiding WF corresponds to a smaller difference between the 
variational Monte Carlo (VMC) and the FN energies. The 
LRDMC energies are computed for a — 0.9 (filled triangles), 
a — 0.5 (open circles) and a — (open squares). The linear 
fits for the DMC and the LRDMC (a = 0.9) data are shown. 
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TABLE I: Comparison of 4s 2 3d n 
ergy (eV) for the scandium atom. 



4s 1 3d n+ excitation en- 



a fully consistent and variational scheme which is much 
more accurate than the standard DMC method. More- 
over, this projection method allows one to deal with 
several length scales through the use of multiple lattice 
spaces, with great reduction of autocorrelation times for 
heavy atoms or complex systems. We believe that this 
framework can have a wide spread of important applica- 
tions ranging from nuclear physics |l8| to the chemistry 
of transition metal compounds. 

This work was partially supported by INFM, by MIUR 
(COFIN 2003), and by Stichting voor Fundamenteel 
Onderzoek der Materie (FOM). We acknowledge G. 
Bachelet and S. De Gironcoli for useful discussions. 



large errors in the DMC total energies, and performs 
the worst for the scandium atom [15j. As before, we 
keep the determinantal part of the WF fixed, and em- 
ploy a 2-body and a 3-body Jastrow factor. The 
determinantal component is an antisymmetrized gemi- 
nal function expanded over a (5s5p5e?) Gaussian-type 
basis in order to cure near degeneracy effects, and op- 
timized in the presence of the 2-body Jastrow factor. 
We employ Dolg's pseudopotential 0] and compute the 
4s 2 3eP — » 4s 1 3gP +1 excitation energy wich is reported 
in Tabic [I] It is apparent that the LA does not only af- 
fect the DMC total energies but also the DMC energy 
differences: the DMC excitation energy computed with 
the 2-body Jastrow factor differs from the experimen- 
tal value by more than three standard deviations. On 
the other hand, the LRDMC FN results, obtained with 
a = 0.5, are much less sensitive to ^g, and are com- 
patible with the experiment even when a simple 2-body 
Jastrow factor is employed. The LRDMC mixed-average 
excitations are also closer to the experimental value than 
the LA ones, probably because the variational treatment 
of non-local pseudopotentials gives a better cancellation 
of errors. Remarkably, we found that a LRDMC simula- 
tion with off-diagonal pseudopotentials is computation- 
ally much more stable and efficient than the standard 
DMC approach since the negative divergences coming 
from the pseudopotentials close to the nodes are con- 
verted to finite hopping terms in the LRDMC scheme. 

We have presented an efficient lattice regularization 
scheme for QMC calculations on realistic electronic sys- 
tems. The main advantage of the LRDMC approach is 
the possibility to work with non-local potentials within 
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